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ABSTRACT 

Over the last two decades the Andromeda Galaxy (M31) has been something of a test-bed for 
methods aimed at obtaining accurate time-domain relative photometry within highly crowded 
fields. Difference imaging methods, originally pioneered towards M31, have evolved into so- 
phisticated methods, such as the Optimal Image Subtraction (OIS) method of Alard & Lupton] 
( |1998| ), that today are most widely used to survey variable stars, transients and microlensing 
events in our own Galaxy. We show that modern difference image (DIA) algorithms such as 
OIS, whilst spectacularly successful towards the Milky Way bulge, may perform badly to- 
wards high surface brightness targets such as the M31 bulge. Poor results can occur in the 
presence of common systematics which add spurious flux contributions to images, such as 
internal reflections, scattered light or fringing. Using data from the Angstrom Project mi- 
crolensing survey of the M31 bulge, we show that very good results are usually obtainable 
by first performing careful photometric alignment prior to using OIS to perform point- spread 
function (PSF) matching. This separation of background matching and PSF matching, a com- 
mon feature of earlier M31 photometry techniques, allows us to take full advantage of the 
powerful PSF matching flexibility offered by OIS towards high surface brightness targets. 
We find that difference images produced this way have noise distributions close to Gaussian, 
showing significant improvement upon results achieved using OIS alone. We show that with 
this correction light-curves of variable stars and transients can be recovered to within '^10 
arcseconds of the M31 nucleus. Our method is simple to implement and is quick enough to 
be incorporated within real-time DIA pipelines. We also demonstrate that OIS is remarkably 
robust even when, as in the case of the central regions of the M3 1 bulge, the sky density of 
variable sources approaches the confusion limit. 
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C3 1 INTRODUCTION 

Difference Image Analysis (DIA) is now used routinely to provide 
very accurate relative photometry of variable stars, transient objects 
and microlensing events in the Milky Way and other nearby galax- 
ies ( |Wozniak|2008| . DIA permits very accurate relative photometry 
even within extremely dense stellar fields where conventional pho- 
tometric methods can fail or suffer from serious bias. 

Most DIA pipelines currently in use derive from the Optimal 
Image Subtraction (OIS) algorithm of [Alard & Lupton|(T998] ) and 
|Alard| ( [2000| . The OIS algorithm has found widest use in provid- 
ing accurate relative photometry of stars in the Milky Way and the 
Magellanic Clouds. Similar schemes that pre-date OIS were orig- 
inally employed to look for microlensing towards the Andromeda 



Galaxy ( [Tomaney & Crotts|[T996l [Ansari et aLl[T997l ). Since the 
stars in the Andromeda Galaxy (M31) are two orders of magnitude 
further away than those typically monitored in our Galaxy, differ- 
ence imaging towards M3 1 throws up several additional challenges 
to the standard technique. 

In this paper we show that towards the bulge of M31, and 
similarly towards other targets where diffuse background surface 
brightness dominates the total flux, DIA pipelines based on the OIS 
algorithm can often yield poor results due to common image sys- 
tematics such as internal reflections, scattered light or fringe ef- 
fects. Systematics, which may appear at a low level 1%) on 
the original exposures, can give rise to large- amplitude differential 
background residuals on difference image frames. Since OIS mini- 
mizes mismatches in both the point spread function (PSF) and the 
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differential background simultaneously, poor background match- 
ing often results in poor PSF matching and therefore substantial 
systematic errors in differential photometry. We propose a straight- 
forward remedy to allow the effects of such systematics to be min- 
imized. 

The structure of the paper is as follows. In Section|2]we briefly 
describe how M3 1 has been used as a test-bed in developing time- 
domain photometry towards crowded stellar fields. The evolution of 
these methods culminated in the Optimal Image Subtraction (OIS) 
algorithm of Alard & Lupton ( 1998 ) and Alard ( 2000 ), which forms 
the basis for most DIA pipelines currently in use. We use images 
obtained by the Angstrom Project ( [Kerins et al.|2006| ) of the bulge 
region of M3 1 to show how the OIS algorithm may not perform op- 
timally in the presence of bright backgrounds. In Section|3]we show 
how difference images with noise levels close to the photon noise 
limit can be recovered by separating the photometric alignment and 
PSF matching stages. In Section|4]we show some example periodic 
variable light-curves from the Angstrom Project dataset, illustrat- 
ing the impact of the correction on their photometric quality and on 
the ability to characterise variable stars at a range of distances from 
the M31 core. We discuss the findings in Section|5] 



2 DIFFERENCE IMAGING TOWARDS THE M31 BULGE 

2.1 Relative photometry in very crowded fields 

The Andromeda Galaxy has provided something of a test-bed for 
the development of algorithms designed to obtain accurate rela- 
tive photometry within crowded fields. At around the same time 
two techniques emerged, difference imaging fTomaney & Crotts] 
^1996 ) and super-pixel photometry ( Ansari et al. 1997 ). Both tech- 
niques deal with the difficulty of obtaining robust relative photom- 
etry across epochs in the presence of variations in seeing and sky 
background. 

|Tomaney & Crotts] ( |199 6) advocated the difference imaging 
approach, in which they convolve a high quality reference image 
with a Gaussian kernel to match the point spread function (PSF) of 
a target image. Prior to the convolution step the smooth M3 1 back- 
ground light was subtracted from both images using a local median 
filter. The PSF matching approach forms the basis of most current 
difference imaging codes, though it assumes the basic functional 
form of the PSF is known a priori. Modern implementations do not 
include a separate background subtraction stage but instead solve 
for the PSF and background simultaneously. 

In the super-pixels method ( |Ansari et al . 1997 ) an image pair is 
also photometrically aligned in order to match the background lev- 
els. This was done first through linear photometric alignment and 
then by local median filtering the images to produce smooth back- 
ground maps. The background map was subtracted from the target 
image and replaced with that of the reference image. In the absence 
of seeing variations, or intrinsic source variability, the distribution 
of image flux deviations (measured within square pixel arrays, or 
super-pixels) with respect to the local background flux should be 
statistically the same between aligned image pairs. By plotting the 
two distributions against one another and determining the best lin- 
ear fit, a simple linear flux correction can be made to the target im- 
age flux to correct for seeing changes. The super-pixels technique 
has the advantage of being empirically calibrated, avoiding strong 
assumptions on the form of the PSF. However in principle it is not 
as sensitive as difference imaging due to the potential loss of signal 
through binning flux into super-pixels. 



Interestingly both of these methods separate out the back- 
ground correction from the seeing correction. Modern difference 
imaging algorithms make these two corrections simultaneously. 
Whilst this is a highly efficient approach, as we shall show, it does 
not always yield optimal results. 

2.2 Optimal image subtraction 

Most difference image packages currently in use employ the Opti- 
mal Image Subtraction (OIS) algorithm present ed by|Alard & Lup-| 
|ton| ( [T998] ) and generahzed further by " Alard] (|2000l). In the OIS 
method a reference image R is convolved to match the PSF of a 
target image T that has been geometrically aligned and re- sampled 
to the same pixel grid as R. By taking into account differences in 
the background between R and T via a smooth 2D background 
model, B, the difference image, D, can be computed from linear 
least-squares minimization: 

D{xi,yif = min{^[(i?(xi, 2/0 K{u, v)) - 

i i 

T{xi,yi) + B{xi,yi)f/a{xi,yif} (1) 

where (x, y) denote the position of image pixel i, K is a. kernel 
function that describes the PSF transformation between R and T, 
(n, v) are coordinates centred on the kernel, is the pixel variance 
and denotes convolution. For computational efficiency |Alard &| 
|Lupton|fl998| ) advocate decomposing K into a linear combination 
of basis functions 

K{u, v) = Y^ ap,iz^i;'^e-("'+"'^/'"' , (2) 

k — l p,q—0 

where Cipq are coefficients. The Gaussian widths ak and integers N 
and Mk control the complexity of the kernel shape. Similarly the 
differential background is modelled as a sum of polynomial basis 
functions: 

r,r-\-s<Mjj 

B{x,y)= ^ arsx'^y% (3) 

r,s—0 

where 

Cirs ^re coefficients and the integer Mb is the degree of poly- 
nomial used to model the differential background. 

There are a number of freely available software implementa- 
tions of the OIS algorithm, the best known of which are ISiaHand 
DIAPlQ Throughout this paper we use ISIS version 2.2 ( |Alard| 
|& Lupton,1998i |Alard|[200 0). However, we stress that the short- 
comings we illustrate are not specific to ISIS but are shared by all 
implementations of OIS that assume simple polynomial forms for 
the differential background as in Equation ([5]). As pointed out in 
the review of |Wozniakl ( |2008| ) the OIS method is actually indepen- 
dent of the choice of basis function. ,Bramich, i^2WS ) has imple- 
mented a version of OIS (DANDIA) in which Equation ([TJ is min- 
imized pixel-by-pixel (i.e. essentially using a (5-function kernel). 
This method solves for the kernel for an assumed constant back- 
ground within a sub region. The global differential background so- 
lution is then found via interpolation over a grid of sub-regions. 
However for complex differential backgrounds of the kind investi- 
gated in this paper the method would potentially require splitting 
the image into a very large number of sub-regions in the absence of 

jhttp : / /www . iap . fr /users /alard/package . html 
^ Ihttp : //users . camk . edu . pl/pych/DIAP?7 
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Figure 1. Left panel: The position of the Angstrom Project robotic field shown on an NOAO image of the Andromeda Galaxy (M31 - image credit: Bill 
Schoening, Vanessa Harvey/REU program/NOAO/AURA/NSF). Right panels: (a) A Sloan i^-band stacked image of the Angstrom Project robotic field ob- 
tained from the Liverpool Telescope at La Palma. This image is used as our reference image for the tests in this paper and covers the inner ~ 4^ x 4^ of the 
M31 bulge, (b) The locations of thousands of variable stars around the M31 core are revealed by combining a stack of 480 difference images. The difference 
images span 5 years and are created from a sequence of more than 3300 individual exposures. The difference images were created using the modified method 
described in this paper. Composite residuals from imperfect subtraction within the inner ~ 10^^ of the core and around some bright foreground stars are also 
visible. 



prior knowledge on the shape of the differential background and its 
behaviour with respect to time. 



2.3 Difference imaging of Angstrom Project data 

We have applied the OIS method to a series of Sloan i^-band im- 
ages of the M31 bulge obtained by the Angstrom Proje ct ([Kerins et 
al.|2006| ) using the Liverpool Telescope (hereafter LT, [Steele et ah 
2004| l on the Canary Islands. The Angstrom Project is undertaking 
long-term monitoring of the bulge of M3 1 in order to detect stellar 
microlensing events and transients. The survey employs a real-time 
difference imaging pipeline capable of issuing alerts of ongoing 
events ( Darnl ey et al.|2007] ). The location of the Angstrom robotic 
telescope field is shown in Figure [T] 

The LT is a 2m robotic telescope housing an optical camera 
(RATCAM) with a 4.'6 x 4.'6 field of view with 0Vl3 pixels. Im- 



ages were obtained over a 5 -year period from 2004-2009. The LT 
typically obtained one or two epochs of data per clear night dur- 
ing the M31 observing season (August-January) with each epoch 
comprising a run of between 7 and 15 short exposures, typically 
ranging from 30 to 200 seconds. The exposures were kept short 
to avoid saturation of the M31 core and bright foreground stars. 
Pre-processing, comprising de-biasing and flat fielding, was car- 
ried out by the LT robotic telescope pipeline. The Angstrom Project 
pipeline then performed preliminary geometric alignment and de- 
fringing, followed by optimal image alignment using Fourier cross 
correlation (for further details see |Darnley et al.|2007| ). For the real- 
time pipeline the images are re-binned to 1024 x 1024 pixels with 
a resolution of 0^.^26. In this paper we also use these re -binned im- 
ages. 

All images were built by aligning and stacking individual ex- 
posures within each observing epoch to form a high signal-to-noise 



© 0000 RAS, MNRAS 000, 000-000 



4 E. Kevins et al. 



epoch frame. The master reference frame was constructed from a 
stack of 30 high quaUty individual exposures obtained over four 
separate epochs from the 2007/8 and 2008/9 observing seasons. 
These epochs were chosen on the basis of accurate telescope point- 
ing and tracking, low airmass, low background levels and good 
seeing. ISIS can itself be used to build the reference image from 
a stack of high quality frames by matching their PSFs and back- 
ground using Equation ([T]). We choose instead to create our ref- 
erence frame by first photometrically aligning the selected images 
(using the method described in Section [3T| l and then median stack- 
ing them using the IRAlQiMCOMBlNE task. Hence our difference 
image results rest only on a single application of Equation ([TJ when 
matching the reference and target images. The same reference im- 
age is used throughout so that variations in our results depend only 
on how we manipulate the target frame. The reference frame is 
shown in panel (a) of Figure [T] 

Figure [2] shows some example difference images created from 
the reference image and a very high quality target image obtained 
under excellent seeing. The results illustrate the challenge of ac- 
curate background matching towards the core of M3 1 where, even 
within the small LT field of view, the dynamic range in surface 
brightness across the image is around two orders of magnitude. In 
Figure |2ja) ISIS uses a second-order polynomial to model the dif- 
ferential background across the whole image area. Whilst the re- 
sulting difference image clearly succeeds in revealing a wealth of 
variable stars (evidenced by the black and white spots, which show 
sources that have dimmed or brightened) it suffers from a signifi- 
cant large-scale background residual. We find that the amplitude of 
the residual is largely insensitive to the choice of polynomial order, 
though the residual pattern across the image does vary according 
to the order used, as one would expect. The underlying differen- 
tial background is therefore clearly not well represented by a 2D 
polynomial. There are a number of possible causes for such an ef- 
fect. Scattered light, internal reflections, dust and fringing effects 
can add localized unwanted flux contributions, whilst flat fielding 
errors may have a multiplicative effect. When observing a target 
like the M31 bulge which has a very bright and spatially varying 
background, the effects of such systematics are unlikely to be ade- 
quately modelled by a polynomial function. 

In fact the amplitude of the background variation over most 
of the image area in Figure |2|a) typically represents only ^ 1% 
of the original image flux, but this is still too large to allow re- 
liable difference image photometry of sources, which themselves 
are varying at a level of a few percent of the background surface 
brightness flux. Since the residual background is epoch dependent 
the failure to effectively remove it means that robust relative pho- 
tometry between epochs becomes virtually impossible. Addition- 
ally, the large- amplitude background residual may blunt the OIS 
algorithm's sensitivity to the optimal PSF transformation because 
the overall summed squared difference flux [D^ in Equation 
can become dominated by the background residual rather than by 
PSF kernel residuals. 

Whilst these effects are important for high surface bright- 
ness regions such as the M3 1 bulge they are far less important for 
Galactic surveys that employ difference imaging. Even towards the 
Galactic bulge the flux observed by optical surveys is dominated 



^ IRAF is distributed by the National Optical Astronomy Observatories, 
which are operated by the Association of Universities for Research in 
Astronomy, Inc., under cooperative agreement with the National Science 
Foundation. 



by resolved or semi-resolved stars and therefore any background 
residuals from systematic problems will be typically at a low level 
even on difference images. However, it is feasible that photometric 
corrections might benefit future near-infrared surveys of the Galac- 
tic bulge, which can target the centre of the Galaxy where the unre- 
solved stellar background flux is considerably higher. Near-infrared 
array technology is also still maturing and so is currently more 
prone to systematic problems than optical CCDs. 

The ISIS software also provides the option of running the OIS 
algorithm independently within sub-regions of an image. In Fig- 
ure [2|b) we show the difference image result where ISIS has split 
the image into 4x4 sub-regions. All other DIA parameters are the 
same as for Figure|2ja). This time ISIS does a better job over much 
of the image area. However, within the inner bulge region, cover- 
ing around half the field area, noticeable background residuals are 
still evident, especially at the square sub-region boundaries, which 
show sharp discontinuity. One might be tempted to try to reduce 
the effect by further increasing the number of sub-regions. How- 
ever, we are fundamentally limited by the size of the image stamps 
within which the convolution kernel is determined, which in turn 
is limited by the scale of the image PSF. For the image in Fig- 
ure [2|b) each sub-region is only around 200 pixels wide, compared 
to a stamp width of 30 pixels. Ideally the PSF and background 
should be evaluated within several independent stamps across the 
sub-region in order to facilitate a good solution for the PSF trans- 
formation. The sub-regions also should be large enough to contain 
resolved or partially resolved stars, or strong features such as dust 
lanes, so that the solution to equation ([T| is not dominated by noise. 

These limitations force us to consider treating the differential 
background matching separately from PSF matching. 



3 A MODIFIED PROCEDURE FOR DIFFERENCE 
IMAGING 

We have seen towards high surface brightness targets such as M31 
systematic effects such as scattered or reflected light can imprint 
complex time- variable signatures upon the differential background 
that the OIS algorithm cannot deal with adequately. The OIS algo- 
rithm performs PSF and background flux matching simultaneously, 
whilst prior to OIS these were treated separately by various M31 
variable photometry pipelines ( Toma ney & Cr otts 1996, Ansari et| 
[aL||1997| ). However, the OIS algorithm provides clear advantages 
over earlier schemes due to its flexibility in modelling the PSF. 
Ideally we would like to combine the best of the old and current 
approaches. 

To this end we choose to separate out the photometric and PSF 
matching stages, as in earlier approaches. Then we run ISIS with 
differential background matching effectively turned off by setting 
the background polynomial to order zero [i.e. setting Mh = in 
Equation |3])]. In this case the least squares minimization of in 
Equation Q should be driven by the quality of the PSF transfor- 
mation kernel K rather than potentially having to trade between 
the quality of PSF and background matching. 

3.1 Photometric alignment 

We begin by determining a gross linear photometric scaling to 
match the image flux. Due to large telescope pointing errors in the 
first two seasons of the Angstrom survey the sky area common to 
all exposures represents less than half the available image area. We 
therefore use the largest available rectangular region common to 
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Figure 2. Upper panels: difference images produced by ISIS (version 2.2) from a pair of high quality reference and target images covering an area of 3.^9 x 3.^7 
around the M31 core. In (a) ISIS is run on photometrically unaligned images allowing a second-order polynomial to model the differential background. Strong 
background residuals remain in the difference image, which do not markedly improve with the use of higher order polynomials. In (b) the photometrically 
unaligned images are split into 16 sub-regions on which ISIS is run individually using the same parameters as in (a). The background matching is better 
further out from the core (apparently at the expense of poorer PSF matching) but remains poor closer to the core where the boundaries between the sub-regions 
are evident. In (c) the images have been photometrically aligned using the method discussed in this paper, prior to running ISIS. In this case a zeroth-order 
polynomial has been imposed for the differential background when running ISIS. All other parameters are the same as in (a) and (b). Lower panels: show a 
1' X 0.^9 zoom of the region shown by the box in (a) for each case (a), (b) and (c), respectively. In (c) note how the removal of the differential background also 
improves the quality of PSF matching as evidenced by the lessened residuals around bright foreground stars. 




Figure 3. Matching the non-linear background residual. A target image (a) is linearly matched to the reference image [shown in Figure [TJa)] via a fit to the 
pixels within the rectangular boundary, a region common to all exposures. The resulting linearly aligned image is then directly subtracted from the reference to 
reveal a residual map (b). The map reveals a large scale non-linear variation in the residual background level. This map is then Gaussian filtered to produce a 
smoothed residual background map (c). The contrast level in (b) and (c) has been stretched by a factor 30 over that in (a) to highlight the background residual, 
which has an amplitude that is typically only about 1% of the target image flux. The smoothed residual map is added back on to the target image to produce 
the final photometrically corrected target frame. The difference image resulting from the corrected target frame is shown in Figure|2jc). 
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X coord, (pixels) X coord, (pixels) 

Figure 5. The reference image and target "structure map" aligned to the same global pixel coordinates, (a) The reference image, showing the locations of 31 
objects (shown circled - mostly foreground stars or features relating to dust lanes) selected by ISIS to determine the convolution kernel. One object is selected 
within each of the square sub-regions which form an overlapping mesh indicated by the dotted lines. In (b) we show a noise-normalised map of residual 
structure in the target image. The structure map is formed by subtracting a smoothed version of the target and then dividing the result by the square-root of the 
target flux. Foreground stars show up clearly on the structure map, as do dust lanes close to the M31 core and in the upper right of the image. The granular 
structure in the image is not noise but partially resolved bright stars and unresolved clumps of stars in M3 1 . The smoothing aperture produces the dark ringing 
effect evident around some of the brighter foreground stars and the M31 core, as well as the border around the image perimeter. 




10000 

Counts 

Figure 4. Flux count histograms of the reference image (points with er- 
rors) and target image (dashed line) after the target has been photometri- 
cally aligned to the reference. The pseudo errors on the reference image 
pixel numbers assume an expected Poisson-like variance in the pixel num- 
bers between aligned images. After alignment the behaviour of the target 
count histogram accurately mirrors that of the reference image even down 
to small scale deviations. 



all images, which covers around 30% of the image area (see Fig- 
ure [3]) in order to photometrically align our images. We determine 
a gross linear scaling and offset between the reference R and target 
T {R = a xT + h) using the LINMATCH task within IRAF. The 
scaling a and offset b are then applied to T to obtain a linearly pho- 



tometrically matched image Tl. After linear scaling we subtract Tl 
from R to obtain a residual image R — Tl. Since the seeing is gen- 
erally different between R and Tl the differential flux from gen- 
uinely variable sources is significantly smoothed out by differences 
in image PSF. A Gaussian filter is then passed over the residual 
image to produce a smoothed residual background image B' . The 
size of the Gaussian is set to a = 15 pixels, which is large enough 
to be insensitive to the presence of variable stars but small enough 
to provide a good local estimate of the differential background. Fi- 
nally a background corrected target image = Tl + B' is ob- 
tained. Figure[3]illustrates how the differential background map B' 
is produced from a target image which has already been linearly 
matched to the reference frame. 

Figure |4] shows the resulting pixel flux histograms after pho- 
tometric matching. Around the M31 bulge the image flux is domi- 
nated almost everywhere by the smooth unresolved surface bright- 
ness distribution, so the photometric calibration is insensitive to 
seeing and we are able to obtain very good photometric alignment 
without the need for PSF corrections. 



The adjusted image T' now has a background that is properly 
matched to R so that ISIS can be run without the inclusion of the 
B term in Equation ([TJ. Note that this matching is effective only 
against systematic effects which are additive (e.g. scattered light). 
One might imagine that a poorly constructed flat field could cause 
a non-uniform multiplicative variation across the image, in which 
case the method described here would not adequately correct the 
flux. 
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Figure 6. (a) Grey-scale "signal-to-noise" map produced by dividing the 
difference image shown in Figure [2jc) by the square root of the combined 
reference and target image counts. Note how the increase in image noise 
towards the M31 core suppresses the number of detectable variations there 
(compare the contrast of sources close to the core with those further out) 
though some are still evident to within ~ 10^^ of the core, (b) A map of all 
the pixels in (a) that exceed 2.5 cr. Isolated random noise pixels are com- 
paratively rare and most belong to PSF-like clusters that are associated with 
genuinely variable objects. These provide the excess seen in the tails of the 
noise histogram shown in Figure^ 



3.2 PSF matching and difference image results 

Before running ISIS on our photometrically aligned data we make 
checks on whether there is enough information on the PSF within 
our images for ISIS to make a reliable PSF kernel determination. 
Since the surface brightness of the M31 bulge is high the number 
of obvious resolved stars is relatively small. 

ISIS determines the PSF kernel transformation not by using 
all the image pixels but by sampling the image within a number 



Error normalized DIA counts 

Figure 7. Histogram of the counts for the image shown in Figure [6ja) to- 
gether with a unit- width Gaussian model. The amplitude of the Gaussian 
is fitted to data points within ±1 cr. The excess from pure Gaussian noise 
evident in the tails of the histogram is predominately associated with gen- 
uinely variable objects as shown in Figure|6] Also shown are the equivalent 
histograms for the difference images shown in Figure|2ja) and (b) where no 
photometric alignment is performed prior to running ISIS. 

of small regions ("stamps") dispersed over the image. The location 
of the stamps can be partly constrained by the user specifying a 
grid of sub-regions. ISIS determines the most suitable stamp loca- 
tion within each sub-region, usually centring the stamp on a bright 
star. The quality of the PSF transformation is therefore partly deter- 
mined by the number of user-defined sub-regions and by the suc- 
cess with which ISIS finds suitable stamps within each sub-region. 

Figure [5] shows the reference image (a) and "structure map" 
of our test target image (b) on which is over-plotted the sub-region 
grid adopted for the tests in this paper. The structure map is pro- 
duced by subtracting a Gaussian smoothed version of the target 
image from the original target frame and then dividing the result 
by the square root of the target pixel flux. This removes the smooth 
background light component of the M3 1 surface brightness, clearly 
revealing foreground Milky Way stars as well as dust lanes, semi- 
resolved bright stars and stellar clumps in M3 1 itself. The structure 
map clearly shows that there is a wealth of structure over the whole 
image area which is potentially useful for PSF determination. 

The grid overlay in Figure [5|a) and (b) is defined with respect 
to a 1536 X 1536 pixel area which defines a global pixel coordinate 
system onto which all images have been aligned. (The large relative 
size of this area is necessitated by large pointing offsets in some of 
the earlier exposures.) The circles in Figure [5|a) and (b) indicate 
the locations of the centre of each of 31 stamps which ISIS selected 
for our test target frame. We set the size of each stamp to 30 x 30 
pixels, so that collectively the stamps comprise about 3% of the 
available image area. More often than not the stamps are centred 
on bright foreground stars, however it is clear from Figure [Sjthat a 
number of the stamps are not centred on any obvious star. We find 
that the steep surface brightness gradient is a strong determinant as 
to whether the ISIS stamp finding algorithm locates a bright star or 
not. 

In order to determine the PSF kernel reliably, the number of 
stamps should be at least 3(1 + M^) (2 + Mfc)/2 ( A lard & Lupton| 
1998, Alard 2000 ), where Mk is the user- specified polynomial or- 
der allowed for the spatial variation of the kernel (see Equation [2]). 
For our data we find Mk = 1 sufficient to characterise spatial PSF 
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kernel variations, so that we require a minimum of 9 good stamps. 
In fact, typically the ISIS stamp finding algorithm results in ^ 15 
stamps centred on genuine bright foreground stars, more than suf- 
ficient to compute a reliable kernel. We have also checked that us- 
ing an improved star finding algorithm, in which all the stamps are 
centred on obvious stars, makes little difference to the overall dif- 
ference image result. In fact, whilst the bright foreground stars are 
preferred for PSF determination, the granular structure of partially 
resolved M31 stars evident in Figure[5|b) means that just about any 
position on the image will provide some useful information for the 
PSF kernel determination. We therefore proceed with the default 
stamp finding algorithm which is contained in the ISIS package. 

The result of running ISIS on a target image that is now pho- 
tometrically aligned to the reference frame is shown in Figure [2jc). 
All ISIS configuration parameters are set the same as for the pre- 
vious tests in Figure [2ja) and (b) except we now enforce a zeroth- 
order polynomial model for the differential background. The re- 
sulting difference image is clearly superior to the previous exam- 
ples. There is now no noticeable background residual anywhere on 
the image. Equally important is the fact that the PSF transforma- 
tion is clearly much improved too. This is evident when examining 
the residuals associated with the few bright foreground stars visi- 
ble within the Angstrom field (compare with the reference image 
shown in Figure [TJ. These foreground stars leave relatively strong 
residuals in the difference images shown in Figure [2ja) and (b) but 
they are much reduced in Figure [2jc) (see also the zoomed inset 
panels). This indicates that ISIS has been able to find a better PSF 
solution as a result of minimizing Equation ([T| for photometrically 
aligned images. 

Figure |6ja) shows a map of the signal to noise (S/N) ratio. 
The signal-to-noise map is defined by 



S/N: 



^/WN{R - T) 
VMT + NR ' 



(4) 



where M is the number of individual exposures comprising the tar- 
get frame stack T, N is the number of individual exposures com- 
prising the reference stack R and g is the CCD gain. The grey- scale 
in Figure |6j a) is calibrated to show pixels with \S/N\ > 10 as ei- 
ther white or black depending upon whether the difference flux is 
positive or negative. The main difference in the appearance of this 
S/N image with the corresponding difference image shown in Fig- 
ure|2jc) is that near the M31 core the S/N image shows much less 
variation, indicating that many of the flux deviations on the differ- 
ence image near the M31 core are just due to the noise arising from 
the very high surface brightness level. Figure |7] shows the result- 
ing pixel histogram of S/N values for the photometrically aligned 
target frame and for the unaligned cases shown in Figure |2| a) and 
(b). These histograms demonstrate a clear improvement both in the 
symmetry and width of the noise distribution as a result of carrying 
out separate photometric alignment. Over-plotted on the noise his- 
togram of the aligned data is a unit-width Gaussian model, where 
the amplitude has been fitted within ±1 a. The histogram shows 
that within this region the noise is well approximated as Gaus- 
sian. Towards the tails of the distribution there is an excess over 
and above Gaussian noise. However, as shown in Figure [6f a) and 
(b), this excess is associated with genuinely varying objects which 
make up a non-negligible fraction of the image pixels. Figures [6] 
and |7] therefore suggest that photometrically aligned images can 
allow high quality difference imaging even to within a few arcsec- 
onds of the M31 bulge. 



3.3 The effect of variables 

The high crowding levels of genuinely variable sources evident in 
Figure [2] potentially poses a problem for the OIS algorithm. The 
minimization in Equation ([T]) would be expected to correspond to 
the optimal difference image only in the limit that, in the absence 
of seeing and background variations, most of the image flux does 
not vary with time. In theory it could be the case that the relative 
phases and amplitudes of variables within a stamp at a given epoch 
could pose a significant source of time variability in the average 
difference flux within the stamp. If so this could also affect the 
residual background model determination. 

In order to suppress the effects of genuine variations on the 
computed kernel ISIS uses sigma clipping both of pixel values 
within stamps and of the distribution of stamps themselves. How- 
ever, for the M3 1 bulge, where a relatively high fraction of pixels 
change due to astrophysical variation, sigma clipping can lead to a 
significant fraction of the image area being excluded, particularly 
around the M31 core. 

To test the potential impact of the high sky density of variables 
we perform a double-pass iterative minimization of Equation ([TJ 
that corrects the target image for the presence of variable sources. 
On the first pass ISIS is run in the normal way. Then, all pixel fluxes 
on the resulting difference image above a specified absolute flux 
threshold are added back on to the target frame to form a modified 
target image in which intrinsically variable sources are largely re- 
moved (i.e. matched to their flux values on the reference image). A 
second pass of ISIS is then performed on this modified target frame. 
The resulting kernel solution for the modified frame is then applied 
to the reference image and the convolved reference is subtracted 
from the original (unmodified) target frame to create the final dif- 
ference image. This method ensures that the final difference image 
is largely insensitive to the presence of variable stars and that as 
much of the image area as possible is used to compute the kernel. 

Reassuringly, we find virtually identical results are recovered 
using either a standard ISIS run or the double-pass method de- 
scribed above, showing that ISIS is capable of finding an optimal 
kernel solution even for high crowding densities of variable stars. 



4 EXAMPLE LIGHT-CURVES 

In Figures [8] and [9] we show two periodic variable star light-curves 
from the Angstrom Project data processed with and without the 
photometric alignment stage described in Section |3.1| Figure |8] 
shows the light-curve of a variable with a period of 501 days lo- 
cated 2'A from the M31 core monitored over five observing sea- 
sons. The right-hand and left-hand columns show the results with 
and without photometric alignment, respectively. Panel (a) shows 
the light-curve itself, panel (b) is the computed periodogram, and 
panel (c) is the result of folding the light-curve with the period cor- 
responding to the peak in the periodogram. 

Looking at panel (a) in the right-hand column of Figure [8] we 
see that the light-curve derived from the photometrically aligned 
images appears much more coherent and less scattered than that 
derived from unaligned data in the left-hand panel. Panel (b) shows 
the computed periodogram for periods from 5 to 5000 days calcu- 
lated using the Generalised Lomb Scargle method of |Zechmeister| 
|& Kiir steT ( 2009 ). This generalised method extends the standard 
Lomb Scargle approach by taking account of errors on the flux and 
and not assuming a zero mean flux, making the periodogram com- 
putation more reliable when there are large time gaps in the data 
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Figure 8. A five-season LT i^-band light-curve of a periodic variable star located 2' A from the M31 core. Results using ISIS applied to photometrically 
unaligned images are shown in the left-hand column whilst those from aligned images are shown in the right-hand column. In each case panel (a) shows the 
light-curve, panel (b) shows the spectral density computed using the Generalised Lomb Scargle method of Zechmeister & Kiirster ( 2009 1 and panel (c) shows 
the result of folding the light-curve in (a) by the period corresponding to the maximum in the spectral density shown in (b). A smoothed light-curve template 
constructed using radial basis function interpolation (see main text) of the folded light-curve is indicated by the smooth curve in (a) and (c). In (c) the folding 
period P, the variable half amplitude A, the asymmetry ratio R and the reduced with respect to the template curve (excluding the 6 worst outlying points) 
are also given. The asymmetry ratio R is the ratio of the rise (minimum up to peak) and fall (peak back to minimum) time-scales, as determined from the 
smooth template. Note how in the unaligned case the systematic scatter has resulted in an incorrect period determination. 



(such as for seasonal surveys like ours). The period corresponding 
to the peak in the spectral density is used to fold the light-curve 
data in panel (c). The folded data in the right-hand column clearly 
exhibits well defined periodic behaviour. 

In order to provide a guide on the performance of the photom- 
etry we compute the reduced with respect to model templates 
of the light-curve behaviour. We assume the underlying behaviour 
of the folded data to be smooth on temporal scales much larger 
than the time resolution of the folded data. This may not always be 
the case so the resulting x^ assessment should be conservative in 
the sense of not being an over-optimistic indicator of the true per- 
formance. We construct a template for each light-curve using the 
SciP]^ implementation of the multiquadric radial basis function 
(RBF) approximation, which is a generalisation to the case of noisy 
data of multiquadric RBF interpolation ( Hard y |197T] ). The RBF ap- 
proximation is controlled by two parameters: a shape parameter {e) 
and a smoothness parameter {S). The shape parameter controls the 
way in which the interpolation function weighting varies with Eu- 
clidean distance (norm) from the data points, whilst the smoothness 



^ http://www.scipy.org/ 



parameter controls the level of smoothing over scatter in the data. 
Setting S — ^ would be equivalent to strict interpolation, where the 
template is forced through every data point (clearly inappropriate 
for noisy data). We fix these parameters throughout at £ = 0.25 
and S — 0.05 which, for adequately sampled light-curves, provide 
a good description of trends in light-curve behaviour without ap- 
pearing to over-fit the data; light-curve "structure" below 1/4 of 
a phase is typically unresolved. For the computation of the multi- 
quadric RBF approximation it is desirable to first rescale the flux 
axis so that both the flux and phase axes have similar numeric scales 
(of order unity). This ensures that the interpolation weighting is 
sensitive to displacements along both axes. 

Looking at the left-hand panels of Figure[8]we see that the sys- 
tematic scatter in the light-curve derived from photometrically un- 
aligned data results in a significantly distorted periodogram which 
peaks at the wrong period. Since for our DIA testing we have not at- 
tempted to exclude poor quality images the light-curves all contain 
a small number of obvious outliers. When computing the reduced 
X^ with respect to the RBF template we therefore exclude the 6 
worst outlying points (representing 1.25% of the full dataset) in or- 
der to obtain a fair measure of how well the template represents the 
bulk of the data. The reduced x^ value reported in panel (c) is quite 
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Figure 9. As for Figurejsjbut this time for a variable star located just 19^^ from the M31 core. In the unaligned case displayed in the left-hand column the 
periodogram fails to identify the correct period and instead the peak spectral density reflects merely the season-to-season variation in average flux. 



acceptable for the aligned data but clearly poor for the unaligned 
folded dataset. 

Figure [9] shows a similar situation for a variable star with a 
period of 294 days located just 19 arcsec from the M31 core, cor- 
responding to a projected distance from the core of just 72 pc for 
an M31 distance of 780 kpc. In the unaligned case the peak in the 
periodogram is completely spurious, corresponding to the seasonal 
variation in average flux due to different parts of the variable phase 
being sampled in different observing seasons. Once again the re- 
duced for the RBF template is quite acceptable for the light- 
curve derived from aligned data but poor for the unaligned dataset. 
The examples of Figures [8] and |9] show that careful photometric 
alignment is essential for the accurate characterisation of variable 
signals. 

It is noticeable that the first season of Angstrom data shows 
good consistency with subsequent seasons for the variable star 
shown in the right-hand column of Figure [s] which is located 2'A 
from the M31 core (about a half of the LT field width away). How- 
ever the light-curve of the variable close to the M31 core shown 
in the right-hand column of Figure |9] is much more erratic during 
the first season than in subsequent seasons. Visual inspection of the 
images shows that the majority of images taken during the LT com- 
missioning season in 2004 show evidence of scattered or reflected 
light close to the core. Improvements in light baffling made to the 
LT in early 2005 substantially improved the quality of later data. 

Figure [To| shows six more examples of detected variable stars 
at angular distances from the M31 core ranging from 223 arcsec 



(near the edge of our reference field) down to just 9 arcsec. They 
illustrate good season-to- season consistency in the corrected pho- 
tometry, as evidenced by the smooth continuity of the folded light- 
curve, and also reasonable characterisation of the noise, as indi- 
cated by the size of the reduced values with respect to the RBF 
template. 

The variable in Figure [T0| located at 9 arcsec from the core 
has a predictably very noisy light-curve and is at the threshold of 
detectability; whilst it displays clear periodic behaviour the struc- 
ture of the variation is not well resolved due to noise. We therefore 
determine that the technique works well at least beyond ^10 arc- 
sec from the M31 core. As Figure [6] shows, there are relatively few 
variable objects that are detectable above the noise so close to the 
M31 core but the examples in Figures [9] and [To| demonstrate that, 
where they are detectable, we are usually able to obtain reliable 
photometry even in this extreme high-background regime. 



5 DISCUSSION 

In this paper we have shown that optimal difference imaging in 
regions of very high background levels, such as the bulge regions 
of galaxies, can be severely compromised due to the presence of 
systematics such as internal reflections, scattered light, or fringe 
effects. In these cases difference images created using the Optimal 
Image Subtraction (OIS) algorithm ( |Alard & Lupton|1998[|Alard| 
|2000| ), which is the most widely used difference image algorithm. 
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Figure 10. Corrected light-curves for a sample of six periodic variable stars over a range of distance from the M31 core. The five-season light-curves are 
shown in the panels labelled (a) together with the smooth template determined from the folded light-curve. The folded data is shown in panel (b) below each 
light-curve. The statistics displayed above the folded data is as for Figures[8]and|9] with the addition of Dm31 which is the angular distance of the object from 
the M31 core in arcseconds. 



may exhibit large amplitude background residuals that make reli- 
able relative photometry difficult if not impossible. Fortunately, we 
have shown that OIS is able to give very good results provided the 
images to be differenced are first photometrically aligned prior to 
difference imaging. 

Using the photometric alignment procedure described in this 
paper we find we can produce difference images of the M3 1 bulge 
that are close to photon noise limited. Not only does it minimize 
or eliminate large amplitude background residuals but it also no- 
ticeably improves the quality of the PSF kernel transformation. We 
have tested the modified pipeline by producing several examples of 



periodic variable light-curves from the Angstrom Project survey of 
the M3 1 bulge. The results allow characterisation of periodic sig- 
nals even to within ~ 10 arcseconds of the M31 nucleus. 

The problem we highlight is specific to targets in which the 
background brightness is high and subject to systematic variations 
which have complex spatial signatures. The OIS method is well 
known to cope admirably for stellar fields within our Galaxy, as 
imaged by current optical surveys, where the image flux is domi- 
nated by resolved or semi-resolved stars rather than by the diffuse 
background light. However, future near-infrared time-domain sur- 
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veys of the Galactic Centre could also conceivably benefit from a 
separation of the photometric and kernel matching stages. 
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